##########################################################################
##########################################################################
######
######  Two functions:  one to create (1-alpha)100% equal-tailed
######                  credible interval and a (1-alpha)100%
######                  HPD credible interval for lambda.
######
##########################################################################
##########################################################################

##### (1-alpha)100% equal-tailed credible interval #####

ETCI = function(y, a = 0, b = 0, alpha){
	n = length(y)
	return(qgamma(c(0.025, 0.975), a + n * mean(y), b + n))
}

ETCI(y, alpha = 0.05)

##### (1-alpha)100% HPD interval #####

HPD.h = function(y, a = 0, b = 0, h = 0.1){
	n = length(y)
	apost = a + n * mean(y)
	bpost = b + n
	mode = (apost - 1) / bpost
	dmode = dgamma(mode, apost, bpost)

	## divide by dmode below to get on scale of 0 to 1 ##
	lint = uniroot(f = function(x){dgamma(x, apost, bpost) / dmode - h}, lower = 0, upper = mode)$root
	uint = uniroot(f = function(x){dgamma(x, apost, bpost) / dmode - h}, lower = mode, upper = 10000)$root

	## calculate coverage, should be around (1-alpha)100% ##
	coverage = pgamma(uint, apost, bpost) - pgamma(lint, apost, bpost)

	return(c(lint, uint, coverage, h))
}

HPD = function(h, y, alpha){
	cov = HPD.h(y, h = h)[3]
	res = (cov - (1 - alpha))^2
	return(res)
}

h.final = optimize(HPD, c(0,1), y = y, alpha = 0.05)$minimum
HPD.h(y, h = h.final)